* Reset settings and initialize log file
launch, path("share/logs_levels")

*-------------------------------------------------------------------------------
* Price and Wasserman (2024), "The Summer Drop in Female Employment"
*
* Description: Plot seasonality in levels vs. logs of different outcomes.
*-------------------------------------------------------------------------------


* Prepare sample and run models
*-------------------------------------------------------------------------------

if "$estimate" != "0" {
	* Load our earnings sample
	gzuse if earn_sample == 1 using "$basepath/data/derived/cps_bms_sample.dta.gz", clear
	rename earnings earn1

	* Alternatively assign positive earnings to salaried workers absent without pay
	gen earn2 = earn1
	replace earn2 = earnweek if emp == 1 & hourly == 0 & absent == 1 & paid_absence == 0 & !missing(earnweek)

	* Alternatively assign positive earnings to salaried teachers absent without pay
	gen earn3 = earn1
	replace earn3 = earnweek if emp == 1 & hourly == 0 & absent == 1 & paid_absence == 0 & school == 1 & teacher == 1 & !missing(earnweek)

	* Report pay status among school-based teachers who absent from work in July
	tab paid_absence if school == 1 & teacher == 1 & emp == 1 & absent == 1 [aw = wtfinl]

	* Compute mean earnings in our earnings sample
	gcollapse (mean) earn1 earn2 earn3 (rawsum) earnwt [pw = earnwt], by(female tm)
	tempfile earnings
	save `earnings'

	* Compute EPOP and mean hours
	gzuse "$basepath/data/derived/cps_bms_sample.dta.gz", clear
	gcollapse (mean) emp (mean) hours (rawsum) wtfinl (first) month tmspline* weeks [pw = wtfinl], by(female tm)

	* Merge in earnings
	merge 1:1 female tm using `earnings', assert(1 3) nogenerate
	tsset female tm

	* Take the log of each outcome
	gen lnemp = ln(emp)
	gen lnhours = ln(hours)
	gen lnearn1 = ln(earn1)
	gen lnearn2 = ln(earn2)
	gen lnearn3 = ln(earn3)

	* Loop over men/women
	foreach f of numlist 0 1 {
		* Specifications
		foreach y in "emp" "hours" "earn1" "earn2" "earn3" {
			if "`y'" == "emp" local wt "wtfinl"
			if "`y'" == "hours" local wt "wtfinl"
			if "`y'" == "earn1" local wt "earnwt"
			if "`y'" == "earn2" local wt "earnwt"
			if "`y'" == "earn3" local wt "earnwt"

			* Run specification in levels
			quietly ivreg2 `y' ib5.month tmspline* weeks if female == `f' [aw = `wt'], bw($bandwidth) robust small
			process_estimates, path("logs_levels") model("f`f'_`y'_lvl")

			* Augment with the average value of the outcome as of May
			quietly sum `y' if female == `f' & month == 5 [aw = `wt']
			quietly estadd scalar may_mean = r(mean)
			estimates save "$basepath/models/logs_levels/f`f'_`y'_lvl.ster", replace

			* Run specification in logs
			quietly ivreg2 ln`y' ib5.month tmspline* weeks if female == `f' [aw = `wt'], bw($bandwidth) robust small
			process_estimates, path("logs_levels") model("f`f'_`y'_log")
		}
	}
}


* Prepare estimates
*-------------------------------------------------------------------------------

* Load estimates into memory
load_estimates, path("logs_levels")

* Prepare coefficient labels for each month
make_coeflabels


* Plot the results
*-------------------------------------------------------------------------------

* Plot estimates
foreach y in "hours" "earn1" {
	foreach l in "lvl" "log" {
		* Prepare options
		if "`y'" == "hours" & "`l'" == "lvl" {
			local title "Average weekly hours"
			local ytitle "Difference relative to May (hours/week)"
			local ymin -4
			local ymax 1
			local ylabel "-4(1)1"
			local rescale
		}
		else if "`y'" == "hours" & "`l'" == "log" {
			local title "Log of average weekly hours"
			local ytitle "100 x log change relative to May"
			local ymin -12
			local ymax 3
			local ylabel "-12(3)3"
			local rescale "rescale(100)"
		}
		else if "`y'" == "earn1" & "`l'" == "lvl" {
			local title "Average weekly earnings"
			local ytitle "Difference relative to May ($/week)"
			local ymin -40
			local ymax 20
			local ylabel "-40(10)20"
			local rescale
		}
		else if "`y'" == "earn1" & "`l'" == "log" {
			local title "Log of average weekly earnings"
			local ytitle "100 x log change relative to May"
			local ymin -6
			local ymax 3
			local ylabel "-6(3)3"
			local rescale "rescale(100)"
		}

		* Prepare background shading
		clear
		set obs 14
		gen rlow = `ymin'
		gen rupp = `ymax'
		gen rval = _n - 0.5

		* Plot estimates
		#delimit ;
		coefplot
			(f1_`y'_`l', offset(-.05) recast(connected) color(black) ciopts(color(black)) label("Women"))
			(f0_`y'_`l', offset(+.05) recast(connected) color(gs6) ciopts(color(gs6)) mfcolor(white) label("Men")),
			coeflabels(`coeflabels')
			title("{bf:`title'}")
			xtitle("")
			xlabel(, alternate)
			ytitle("`ytitle'")
			yscale(range(`ymin' `ymax'))
			ylabel(`ylabel')
			addplot(
				rarea rupp rlow rval if inrange(rval, 0.5, 5.5), color($ltgs) plotregion(margin(t=0 b=0)) below ||
				scatteri 0 0.5 0 13.5, recast(line) color(black) below)
			vertical
			`rescale'
			legend(rows(1) size(*.8) label(2 "Women") label(5 "Men"));
		#delimit cr

		tempfile g`l'
		graph save "`g`l''"
	}

	grc1leg "`glvl'" "`glog'", imargin(vsmall)
	nicepdf "$basepath/output/logs_levels_`y'.pdf", panels(2) indirect replace
}

* Report the estimates
esttab f1_emp_lvl f1_emp_log f0_emp_lvl f0_emp_log, transform(@*100 @*100, pattern(1 1 1 1)) compress
esttab f1_hours_lvl f1_hours_log f0_hours_lvl f0_hours_log, transform(@*100 @*100, pattern(0 1 0 1)) compress
esttab f1_earn1_lvl f1_earn1_log f0_earn1_lvl f0_earn1_log, transform(@*100 @*100, pattern(0 1 0 1)) compress
esttab f1_earn2_lvl f1_earn2_log f0_earn2_lvl f0_earn2_log, transform(@*100 @*100, pattern(0 1 0 1)) compress
esttab f1_earn3_lvl f1_earn3_log f0_earn3_lvl f0_earn3_log, transform(@*100 @*100, pattern(0 1 0 1)) compress


* Report absolute and percent summer changes in earnings by sex
*-------------------------------------------------------------------------------

* Loop over outcomes
foreach y in "emp" "hours" "earn1" "earn2" "earn3" {
	if "`y'" == "emp" {
		local ylbl "EPOP"
		local yscale = 100
	}
	else if "`y'" == "hours" {
		local ylbl "hours"
		local yscale = 1
	}
	else if "`y'" == "earn1" {
		local ylbl "earnings"
		local yscale = 1
	}
	else if "`y'" == "earn2" {
		local ylbl "earnings (adjust all absent without pay)"
		local yscale = 1
	}
	else if "`y'" == "earn3" {
		local ylbl "earnings (adjust teachers absent without pay)"
		local yscale = 1
	}

	* Loop over gender
	foreach f of numlist 0 1 {
		if `f' == 0 local flbl "Men:  "
		if `f' == 1 local flbl "Women:"

		* May-to-July changes: absolute, percent change from level spec, percent change from log spec
		quietly estimates restore f`f'_`y'_lvl
		display "`flbl' May to July abs. change in `ylbl': " %5.2f `=`yscale' * _b[coef2]'
		display "`flbl' May to July pct. change in `ylbl': " %5.2f `=100 * _b[coef2]/e(may_mean)'
		quietly estimates restore f`f'_`y'_log
		display "`flbl' May to July log% change in `ylbl': " %5.2f `=100 * (exp(_b[coef2]) - 1)' _n

		* Average change in earnings over the summer
		if inlist("`y'", "earn1", "earn2", "earn3") {
			quietly estimates restore f`f'_`y'_lvl
			local avg_change = 1/3 * (_b[coef1] + _b[coef2] + _b[coef3])
			display "`flbl' Avg. summer abs. drop in `ylbl': " %5.2f `avg_change'
			display "`flbl' Avg. summer pct. drop in `ylbl': " %5.2f 100 * `avg_change'/e(may_mean) _n
		}
	}
}

* Close the log file
unlaunch
